Redox signaling by glutathione peroxidase 2 links vascular modulation to metabolic plasticity of breast cancer

Significance Redox regulation of breast cancer underlies malignant progression. Loss of the antioxidant glutathione peroxidase 2 in breast cancer cells increases reactive oxygen species, thereby activating hypoxia inducible factor-α (HIF1α) signaling. This in turn causes vascular malfunction, resulting in hypoxia and metabolic heterogeneity. HIF1α suppresses oxidative phosphorylation and stimulates glycolysis (the Warburg effect) in most of the tumor, except for one cancer subpopulation, which was capable of using both metabolic modalities. Hence, adopting a hybrid metabolic state may allow tumor cells to survive under aerobic or hypoxic conditions, a vulnerability that may be exploited for therapeutic targeting by either metabolic or redox-based strategies.

(D) Bar graph showing log2[FC] values for the differentially expressed genes-associated with EMT/stemness in cluster 5 comparing GPx2 KD tumor cells (n=520) and control tumor cells (720); *p < 0.05；****p < 0.0001.    Single cell RNA sequencing of the GPx2 KD tumor (g) revealed 7 epithelial tumor cell clusters (h), which, with the exception of cluster 5 (i), exhibited aerobic glycolysis or Warburg effect (j). By contrast, cluster 5 was endowed with the ability to use both OXPHOS and glycolysis, implying metabolic plasticity (k). The latter was demonstrated by the dual expression of p-AMPK (TRITC) and GLUT1 (FITC) in cluster 5-like cells, as markers of OXPHOS and glycolysis respectively (l). By contrast, clusters 0, 1, 2, 3, 4, 6, which exhibit the Warburg effect, are predominantly GLUT1 (FITC) positive (l). These data suggest that cluster 5 may drive tumor cell adaptation and survival, underlying an aggressive phenotype leading to tumor progression.

Cell lines
Primary mammary tumor cell lines were prepared as follows. Tumors were excised, washed in PBS and diced with a razor blade into smaller fragments. Each gram of tumor was incubated with 10 mL of digestion media (1 mg/mL Collagenase type I, 100 unit/ml Hyaluronidase, 100 units/ml Penicillin/ streptomycin, 2 mg/ml of BSA in Medium 199). The suspension was incubated at 37 ºC for 3 hrs with occasional mixing. Digested material was spun at 1000 rpm for 5 min and pelleted cells were plated overnight at high density (1X107 cells per 10 cm dish) in DME/10%FBS supplemented with 10 μg/ml insulin and 20 μg/ml EGF. Fibroblasts were depleted from epithelial tumor cells by limiting trypsinization. Epithelial cells were grown in culture for several weeks till they reached crisis and adapted to in vitro culture conditions. The human MDA-MB-361 and MDA-MB-231 breast cancer cell lines were cultured in DMEM medium (Gibco) supplemented with 10% FBS and 1% Pen-Strep (Gibco).

Animal studies
Animal protocols used for this study were reviewed and approved by the Institute for Animal Studies at

Tumor growth kinetics
Tumor growth curve was determined by measuring tumor size. The tumor size was measured twice a week using a caliper. Tumor volume was determined upon the formula: tumor volume = shorter diameter2 X longer diameter/2. Mice (five per group for PyMT1 vs. PyMT1/GPx2 KD; four for PyMT2 vs five for PyMT2/GPx2-OE; three for MDA-MD-361 vs three for MDA-MB-361/GPx2sh1) were sacrificed at end point.
Data are displayed as mean tumor volume -/+ SEM. Statistical analyses were performed comparing individual time points by unpaired t-test and significant differences were established as p-value < 0.05.

Lung metastasis
Mice bearing control and GPx2 KD or GPx2 OE mammary tumors were sacrificed at end point of 1-2 cm diameter, as allowed by our animal protocol. Lungs were inflated by tracheal cannulation with injection of 1-2 ml of 10% neutral buffered formalin. Formalin fixed lungs were paraffin embedded and blocks sectioned on a tissue microtome (Leica Microsystems) at 5 μm. Lungs were serial sectioned through the tissue and sets of 5 serial sections, at 300 um intervals. Analysis was performed on whole sections after it was determined by inspection that metastases seeded in random locations. Data are displayed as mean foci number ± SEM.
Statistical analyses were performed using unpaired t-test and significance determined at p < 0.05.

GPx2 lentivirus overexpression
The open reading frame of mouse GPx2 (GenScript NM_030677.2) was subcloned by PCR into Biosystems. Murine GAPDH was used as internal reference gene and relative mRNA levels of candidate genes were determined using the 2(-ΔΔCT) method. Three independent experiments were performed.
Statistical significance was calculated using unpaired t-test, p < 0.05.

Immunofluorescence
Cells seeded on coverslips were fixed in 3.7% paraformaldehyde in PBS, permeabilized in 0.01% Triton X-100, washed in PBS and blocked in 2% BSA in PBS for 1hr. Cells were incubated with primary antibody in PBS/0.5% BSA for 2 hrs, washed and incubated with 1:2000 Alexa Fluor secondary antibody for 1 hr. Cells were washed in PBS/0.5% BSA and counterstained with DAPI. For tissue immunofluorescence, formalinfixed/paraffin-embedded tumor tissues were sectioned in 5 μm thickness, deparaffinized in xylene, and rehydrated in a series of 100% ethanol, 95% ethanol, and distilled water. Antigen retrieval was performed 1X antigen retriever solution (Sigma, pH 6.0). tissues were incubated with primary antibody in 5% donkey serum, 2% BSA, 0.5%TX-100 in TBS, followed by incubation of Alexa Fluor secondary antibody for 1.5 hr at room temperature, washed and counterstained with DAPI to visualize nuclei.

ROS Measurements
ROS detection reagent, 2′,7′-dichlorofluorescein (DCF) (Invitrogen, C6827) was used to determine intracellular ROS levels. DCF was reconstituted using anhydrous dimethylsulfoxide (DMSO, Sigma) with optimal working concentration at 10 μM in DMEM/2% serum. Media with DCF were added to cells in 96 wells and incubated for 1 hr in the dark at a 37 ºC, 5% CO2 incubator. After incubation, media with DCF was removed and cells were washed twice with PBS. DCF fluorescence intensity relative to background fluorescence in wells that do not contain cells was determined using a fluorescent microplate reader using excitation wavelength at 495nm and emission wavelength at 520nm. The results were normalized to cells/well. Data are displayed as Mean H2DCFDA fluorescence intensity ± SEM.

In vitro and in vivo cell proliferation
Tumor cell proliferation in vitro was performed using the Premix WST-1 (Takara, Cat. MK400). Cells were seeded in 96 wells, flat bottom wells, in a 100 μl/well culture medium for the indicated time points; 10 μl/well Premix WST-1 was added to each well for 1 hr at tissue culture incubator. WST-1 is tetrazolium salt that is cleaved to formazan dye by succinate-tetrazolium reductase in viable cells, and measured at wavelength 420~480 nm. Tumor cell proliferation in vivo was measured by immunostaining of tumor sections using anti-Ki67 antibody (Abcam, AB15580). Two to three sections per tumor was stained, and 5 random fields per section were taken by using Zeiss fluorescence microscopy for counting Ki67 positive cells. Data are displayed as mean Ki67 positive cells ± SEM.

Transwell Matrigel Invasion
Cell invasion assays were performed in Corning Biocoat Matrigel Invasion chamber according to the manufacturer's instructions. Cell monolayers were made into single cell and plated at 1X105 cells/0.5 ml DME, 0.5% FBS per well into the upper compartment of the Boyden chamber for 18 hrs. Media with 10% serum were used in the lower chamber. Cells on top the filters were removed by cotton swabs, and cells penetrating the filters were stained with 0.5% crystal violet/1% formalin/20% methanol. Invaded cells were photographed and quantified by imageJ. The average number of invaded cells bound per microscopic field over five fields per assay was measured in triplicate experiments. Statistical significance was determined by one-tailed t-test.

Immunohistochemistry
Formalin-fixed/paraffin embedded tissues were sectioned at 5μm, deparaffinized in xylene, and

Vessel Perfusion
Tumor vessel perfusion assay was analyzed in three control and three mice bearing two mammary tumors on bilateral sites resulting in 6 tumors per group. Mice were injected intravenously with 50 μg DyLight® 594 labeled tomato lectin (Vector Laboratories, Cat# DL-1177) in saline. 10 minutes later, mice were euthanized and perfused with 1% PFA. Six mammary tumors were harvested from either control or GPx2 KD tumor bearing mice, and frozen in optimum cutting temperature compound. Frozen sections were stained with anti-endomucin antibody (Santa cruz, sc-65495) that labels all vessels (functional and non-functional).
Secondary antibody (Alexa Fluor 488 donkey anti-rat) was used to visualize endomucin labeled vessels.
Perfusion ratio was calculated by dividing lectin-positive vessels by endomucin positive vessels using imageJ.
Unpaired t-test was used to determine significance; p<0.05.

Vessel Maturation
Vessel maturation was analyzed by co-immunostaining of 6 control and 6 GPx2 KD mammary tumors with anti-endomucin antibody which labels vessels together withanti-laminin which stains basement membranes, or with anti-PDGFRb which stains pericytes. Basement membrane or pericyte deficient vessels (immature) were mostly labeled with endomucin whereas mature vessels were labeled with endomucin and laminin. The maturation fraction of tumor vessels was quantified by ImageJ, and the maturation ratio was calculated by dividing laminin or positive vessels by entire endomucin labeled vessels.

Oxygen consumption in vitro
Cells were suspended in normal growth medium at a concentration of 1 X10 4 cells/ml, and 100 μl of cells was seeded on Seahorse 96-well plates 48 hrs prior to performing the assay. Oxygen consumption rate (OCR) assay was carried out in a XF96 Seahorse Analyzer (Seahorse Bioscience, Billerica, MA, USA). On the day of the assay, medium was replaced with prewarmed (37 °C) 180 μl Seahorse medium (Seahorse XF medium with 2 mM glutamine, 10 mM glucose, 1 mM pyruvate) and the plate was incubated for 1 hour in a 37 °C non-CO2 incubator. The mitostress test was performed with 1 μM oligomycin, 1 μM FCCP, and 0.5 μM rotenone/antimycin. The final OCR values were normalized to cells/well using the CyQUANT Cell Proliferation Assay Kit (Invitrogen, C7026).

Oxygen consumption in vivo
Tissue bioenergetics was determined using an XFe24 Seahorse respirometer. Briefly, tumor tissue was collected rapidly after sacrifice, cut into small pieces (8-12 mg) and quickly transferred to individual wells of an XFe24 plate. Individual pieces were stabilized from excessive movement by islet capture screens (Seahorse Bioscience), and 650μL Krebs-Henseleit buffer (KHB) (111 mM NaCl, 4.7 mM KCl, 2 mM MgSO4, 1.2 mM Na2HPO4, 0.5 mM carnitine, 2.5 mM glucose and 10 mM sodium pyruvate) were added to each well. Digitonin (0.5 μg/mL, Sigma, D6628) was added to enhance plasma membrane permeability. Basal oxygen consumption rates (OCR) were determined at 37°C according to the following plan: Basal readings recorded every 2 min for 5 readings, followed by exposure to digitonin. Subsequent readings were recorded after 2 min mixing and 2 min rest. Basal OCR values were normalized to individual tissue weights.

Glycolysis stress test in vitro
Cells were suspended in normal growth medium at a concentration of 1 x10 4 cells/ml, and 100 μl of cells was seeded on Seahorse 96-well plates 48 hrs prior to performing the assay. Glycolysis stress assay was

Quality Control and Normalization
The unfiltered gene-barcode matrix for each sample was imported into R 3.6.1 and converted into SingleCellExperiment (SCE) objects. Any barcodes with a total UMI count greater than 200 were classified as cells. Seurat 3.1.1 was used for further quality control and preprocessing of the exported SCE objects. Cells with a minimum cut-off of 200 and a maximum cut-off of 4500 total RNA count were filtered. In addition, cells with a percentage of total reads that aligned to mitochondrial genome greater than 10% were removed (inferred as cells undergoing stress and cell death). After quality control, a total of 7621 cells from PyMT1 control tumor and 7768 cells from PyMT1/GPx2 KD tumor were used for normalization and downstream analyses. Cells passing QC were normalized using the SCTransform pipeline in Seurat 3.1.1.

Integration and Clustering
After normalization, the data obtained from PyMT1 control and PyMT1/GPx2 KD tumor underwent comprehensive integration using the integration pipeline in Seurat with the genes commonly retained between control and GPx2 KD tumor serving as integration anchors (4). Principal component analysis was performed to enable unsupervised clustering on the integrated data from PyMT1 control and PyMT1/GPx2 KD tumor, and clustered cells were projected in two dimension using Uniform Manifold Approximation and Projection (UMAP) for cluster visualization (5). Louvain community detection algorithm was further applied to identify the cluster partitions, with a resolution parameter of 0.4 that was identified by building a clustering tree on the datasets to optimize the number of clusters.

Cluster Comparison and Visualization
To compare cluster specific genes under two conditions, we used feature plot and violin plot functions to highlight expression patterns of marker genes of interest. Dot plots were used to visualize the average expression pattern of percentage and intensity of expression of genes of interest cells across all cell clusters.

Ingenuity Pathway Analysis
To investigate how GPx2 knockdown in PyMT1 tumor affects global gene expression, all genes regardless of subpopulations and differential expression status were ranked in decreasing order upon their average difference from the overall differential expression analysis without measuring clusters or cell state specific effect. For bulk analysis, clusters 7 and 8 were removed to get the overall differentially expressed genes. Genes with an adjusted p-value less than 0.05 were classified as differentially expressed. The dataset was uploaded to Ingenuity Pathways Analysis ( IPA) www.ingenuity.com), for core analysis, and a p-value of 0.05 was used as a cut-off to determine statistically significant enrichment of a pathway or annotated gene grouping present in the Ingenuity Knowledge base. A Z-score greater than +2 or less than -2 was used as a cut-off to predict a pathway activation or inhibition, respectively.
To investigate the upstream and downstream effects underlying functional activation or inhibition of OXPHOS or glycolysis pathways, overlay of the differentially expressed genes by GPx2 KD onto genes governing OXPHOS or glycolysis pathways in Ingenuity Knowledge Base allows IPA to determine functional activation or inhibition of oxygen consumption or glucose metabolism.
To investigate gene expression change by cluster comparison between PyMT1 control and PyMT1/GPx2 KD tumor, a joint variable of condition * cell-type was used to calculate cluster-type specific differential expression. Similarly, all differentially expressed genes belonging to specific cluster regardless of differential expression status were ranked in decreasing order based on their average difference. The data was then uploaded to IPA for core analysis.

TCGA Data Mining
The Tumor Cancer Genome Atlas (TCGA) breast cancer RNA-seq data was downloaded from the https://xenabrowser.net. Differences in GPx2 mRNA expression between tumor and matched normal breast tissue was established by paired t-test. Differences in expression of GPx2 mRNA between the four breast cancer subtypes (https://www.nature.com/articles/nature11412) were determined by the Kurksal Wallis test EdgeR (version 3.10) was used to generate a gene rank list of BC subtypes with low and high GPx2 mRNA expression. BC subtypes were divided into quartiles with quartile 1 representing 25% of tumors with the lowest GPx2 mRNA expression and quartile 4 with 25% of tumors expressing the highest levels of GPx2 mRNA.

Kaplan-Meier Analysis
Kaplan-Meier analysis comparing the expression levels of GPx1, GPx2, GPx3, GPx4 mRNA across all four breast cancer molecular subtypes with patient survival duration was performed using a breast cancer data set of 1809 patients using (http://kmplot.com/analysis/index.php?p=service&start=1).

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was used to determine the gene sets enriched in breast cancer subtypes with low expression of GPx2 (https://www.ncbi.nlm.nih.gov/pubmed/16199517). Gene sets with a false discovery rate less than 0.15 were considered significant

Single Cell Sequencing data availability
The single cell RNA sequencing data reported in this manuscript are available at https://www.ncbi.nlm.nih.gov/geo with access number GSE152368. Coding analyses for single cell RNA sequencing data are available in Github with the URL accession link: https://github.com/Malindrie/Breastcancer-scRNA-seq-analysis

Statistical Analysis
Results represent the Mean ± SEM (standard error of the Mean) for indicated experiments with 2-3 independent biological replicas. The statistical methods used are described in the figure legends. Most data were analyzed by the unpaired two-tailed Student's t test to compare two groups with significance set at the p value of < 0.05. Statistical analysis was performed using the Prism 8 (GraphPad) software.